package ru.vsu.amm.calculus;

import java.util.Collections;
import java.util.Comparator;
import java.util.List;

/**
 * <p>
 * СЛАУ с обобщённой матрицей жёсткости параллелепипеда
 * и правой частью, полученной на основании заданных перемещений в узлах.
 * 
 * @author Илья Иванов
 * @since 02 декабря 2012
 *
 */
public class GeneralizedSystem {
	
	/** Размерность матрицы системы */
	private int N;
	
	/** Модуль Юнга */
	private double E;
	
	/** Коэффициент Пуассона */
	private double nu;
	
	/** Шаг сетки по оси x */
	private double hx;
	
	/** Шаг сетки по оси y */
	private double hy;
	
	/** Шаг сетки по оси z */
	private double hz;
	
	/**
	 * Координаты узлов сетки в терминах смещений
	 * { nx, ny, nz} - номер узла соответствует индексу набора смещений в массиве
	 */
	private int[][] crd;
	
	/**
	 * Таблица смежности сетки
	 * { i, j, m, p} - правоориентированная нумерация тетраэдра 
	 */
	private int[][] adj;
	
	/** 
	 * Заданные перемещения в узлах
	 *   int         double double double
	 * { номер узла, u,     v,     w     } - элемент списка
	 */
	private List<Object[]> mvt;
	
	/** Компаратор элементов таблицы заданных перемещений <code>mvt</code> */
	private static final Comparator<Object[]> MOVEMENTS_COMPARATOR = new Comparator<Object[]>() {
		public int compare(Object[] o1, Object[] o2) {
			return ((Integer) o1[0]).compareTo((Integer) o2[0]);
		}
	};
	
	/** Ключ для поиска узла в таблице перемещений */
	private static final Object[] SEARCH_KEY = new Object[1];
	
	/** Ширина ленты */
	private int L;
	
	/** Объём частичного параллелепипеда */
	private double V;
	
	/** Потсоянный коэффициент перед матрицами жёсткости 3 x 3 */
	private double mult;	
	
	/** Постоянная матрица, характеризующая жёсткость материала */
	private double[][] D = new double[6][6];
	
	/** Коэффициенты функции форм */
	private double[] coeff = new double[3];
	
	/** Матрицы перемещений для вершин */
	private double[][] Bi = new double[6][3];
	private double[][] Bj = new double[6][3];
	private double[][] Bm = new double[6][3];	
	private double[][] Bp = new double[6][3];
	
	/** Текущие вычисленные матрицы Bk */
	private Object[] B = new Object[4];
	
	/** Матрица жёсткости ребра rs */
	private double[][] Krs = new double[3][3];
	
	/** Перемещение в узле */
	private double[] kMvt = new double[3];
	
	/** 
	 * Верхняя часть ленты обобщённой матрицы
	 * 0  1 ... L - 1
	 * |  |     |
	 * |  |     |
	 * |  |     |
	 * |  |   N - L + 1
	 * |  |
	 * |  |
	 * | N - 1
	 * N 
	 */
	private double[][] bnd;
	
	/** Правая часть */
	private double[] f;
	
	private GeneralizedSystem() {}
	
	/**
	 * Формирует СЛАУ с обобщённой матрицей жёсткости для параллелепипеда
	 * и правой частью, оперделяемой заданными перемещениями в узлах.
	 * @param E Модуль Юнга (параметр материала), E >= 0
	 * @param nu Коэффициент Пуассона (параметр материала), nu <> 1
	 * @param hx Величина шага сетки по оси x
	 * @param hy Величина шага сетки по оси y
	 * @param hz Величина шага сетки по оси z
	 * @param crd Таблица координат узлов сетки в терминах смещений
	 * @param adj Таблица смежности сетки
	 * @param mvt Граничные условия дифференциальной задачи в узлах сетки
	 * @param L Ширина ленты будущей матрицы (максимальный модуль разности между смежными элементами плюс единица)
	 * @return Обобщённая СЛАУ
	 */
	public static GeneralizedSystem formGeneralizedSystem(
		double E,
		double nu,
		double hx, 
		double hy,
		double hz,
		int[][] crd,
		int[][] adj,
		List<Object[]> mvt,
		int L) {
		
		if (L < 1)
			throw new RuntimeException("Ширина ленты должна быть положительной");
		
		GeneralizedSystem system = new GeneralizedSystem();
		
		system.N = crd.length * 3; // в каждом узле три неизвестных перемещения (u, v, w)
		system.E = E;
		system.nu = nu;
		system.hx = hx;
		system.hy = hy;
		system.hz = hz;
		system.crd = crd;
		system.adj = adj;
		system.mvt = mvt;
		system.L = 3 * L;
		system.V = hx * hy * hz;
		
		system.bnd = new double[system.L][];
		for (int n = system.N, i = 0; i < system.L; n--, i++)
			system.bnd[i] = new double[n];
		
		system.f = new double[system.N];
		
		system.calculate();
		system.clean();
		
		return system;
	}
	
	/**
	 * @return Верхяя часть ленты обобщённой матрицы жёсткости.
	 */
	public double[][] getBand() {
		return bnd;
	}
	
	/**
	 * @return Правая часть обобщённой СЛАУ.
	 */
	public double[] getF() {
		return f;
	}
	
	/**
	 * Основные расчёты.
	 */
	private void calculate() {
		calcMult();
		calcMatrixD();
		
		// Сортируем заданные перемещения в порядке возрастания номеров узлов
		Collections.sort(mvt, MOVEMENTS_COMPARATOR);
		
		// Обходим все конечные элементы - тетраэдры. Делаем разборос эелментов матрицы
		// для каждого тетраэдра в обобщённую матрицу.
		for (int[] tetr : adj) {
			if (tetr[0] >= crd.length || tetr[1] >= crd.length || tetr[2] >= crd.length || tetr[3] >= crd.length)
				throw new RuntimeException("Таблица смежности задана некорректно");
			
			processTetrahedron(tetr);
		}
	}
	
	/**
	 * Обрабатывает один конечный элемент - тетраэдр.
	 * @param tetr Вершины тетраэдра в правильной нумерации
	 */
	private void processTetrahedron(int[] tetr) {
		for (int k = 0; k < B.length; k++)
			B[k] = null;
		
		for (int a = 0; a < tetr.length; a++)
			for (int b = 0; b < tetr.length; b++) {
				int r = tetr[a];
				int s = tetr[b];
				
				double[] rMvt = getMovement(r);
				if (rMvt != null) {
					if (r == s)
						throwOnes(r, rMvt);
					// В противном случае в матрицу и правую часть идут нули - ничего не делаем
					
					continue;
				}
				
				double[] sMvt = getMovement(s);
				if (sMvt != null || r <= s) {
					double[][] Br = calcMatrixBk(tetr[0], tetr[1], tetr[2], tetr[3], a);
					double[][] Bs = calcMatrixBk(tetr[0], tetr[1], tetr[2], tetr[3], b);
					double[][] Krs = calcMatrixKrs(Br, Bs);
					if (sMvt != null)
						throwToRigth(r, Krs, sMvt);
					else
						throwTopTriangle(r, s, Krs);
				}
			}
	}
	
	/**
	 * Вычисляет постоянный коэффициент перед матрицами жёсткости 3 x 3.
	 */
	private void calcMult() {
		mult = E * (1 - nu) / (6 * V * (1 + nu) * (1 - 2 * nu));
	}
	
	/**
	 * Вычисляет постоянную матрицу, характеризующую параметры материала.
	 */
	private void calcMatrixD() {
		int i = 0;
		for (; i < 3; i++)
			D[i][i] = 1;
		
		double d = (1 - 2 * nu) / (2 * (1 - nu));
		for (; i < 6; i++)
			D[i][i] = d;
		
		d = nu / (1 - nu);
		for (i = 0; i < 2; i++)
			for (int j = i + 1; j < 3; j++)
				D[i][j] = D[j][i] = d;
	}
	
	/**
	 * Внимание: тройка векторов ij, im, ip правая!
	 * @param i Вершина тетраэдра 0
	 * @param j Вершина тетраэдра 1
	 * @param m Вершина тетраэдра 2
	 * @param p Вершина тетраэдра 3
	 * @param k Номер вершины тетраэдра, для которой строится функция формы (матрица B)
	 * @return Матрица B для вершины k
	 */	
	private double[][] calcMatrixBk(int i, int j, int m, int p, int k) {
		if (B[k] != null)
			return (double[][]) B[k];
		
		double[][] Bk;
		if (k == 0) {
			Bk = Bi;
			calcFormCoeff(i, j, m, p, coeff);
		} else if (k == 1) {
			Bk = Bj;
			calcFormCoeff(j, m, i, p, coeff);
		} else if (k == 2) {
			Bk = Bm;
			calcFormCoeff(m, i, j, p, coeff);
		} else {
			Bk = Bp;
			calcFormCoeff(p, i, m, j, coeff);
		}
		B[k] = Bk;
		
		zero(Bk);
		Bk[0][0] = coeff[0];
		Bk[1][1] = coeff[1];
		Bk[2][2] = coeff[2];
		Bk[3][0] = coeff[1];
		Bk[3][1] = coeff[0];
		Bk[4][1] = coeff[2];
		Bk[4][2] = coeff[1];
		Bk[5][0] = coeff[2];
		Bk[5][2] = coeff[0];
		
		return Bk;
	}
	
	/**
	 * Внимание: тройка векторов ij, im, ip правая!
	 * @param i Номер вершины тетраэдра
	 * @param j Номер вершины тетраэдра
	 * @param m Номер вершины тетраэдра
	 * @param p Номер вершины тетраэдра
	 * @param coeff Коэффициенты b, c, d функции формы Ni для тетраэлра ijmp.
	 */
	private void calcFormCoeff(int i, int j, int m, int p, double[] coeff) {
		// bi = (zmj * ypj - ymj * zpj) 
		coeff[0] = ((crd[m][2] - crd[j][2]) * (crd[p][1] - crd[j][1]) - (crd[m][1] - crd[j][1]) * (crd[p][2] - crd[j][2])) * hy * hz;
		// ci = (xmj * zpj - zmj * xpj)
		coeff[1] = ((crd[m][0] - crd[j][0]) * (crd[p][2] - crd[j][2]) - (crd[m][2] - crd[j][2]) * (crd[p][0] - crd[j][0])) * hx * hz;
		// di = (ymj * xpj - xmj * ypj)
		coeff[2] = ((crd[m][1] - crd[j][1]) * (crd[p][0] - crd[j][0]) - (crd[m][0] - crd[j][0]) * (crd[p][1] - crd[j][1])) * hx * hy;
	}
	
	/**
	 * 
	 * @param Br Матрица Br 3 x 3
	 * @param Bs Матрица Bs 3 x 3
	 */
	private double[][] calcMatrixKrs(double[][] Br, double[][] Bs) {
		zero(Krs);
		
		for (int p = 0; p < 3; p++)
			for (int q = 0; q < 3; q++) {
				for (int l = 0; l < 6; l++) {
					double d = 0;
					for (int k = 0; k < 6; k++)
						d += Br[k][p] * D[k][l];
					Krs[p][q] += Bs[l][q] * d;
				}
				
				Krs[p][q] *= mult;
			}
		
		return Krs;
	}
	
	/**
	 * 
	 * @param i Номер узла
	 * @return Перемещение в узле i, если оно задано, или <code>null</code> в противном случае.
	 */
	private double[] getMovement(int i) {
		SEARCH_KEY[0] = Integer.valueOf(i);
		int found = Collections.binarySearch(mvt, SEARCH_KEY, MOVEMENTS_COMPARATOR);
		if (found < 0)
			return null;
		
		Object[] kMvt = mvt.get(found);
		for (i = 0; i < 3; i++)
			this.kMvt[i] = (double) kMvt[i + 1];
		
		return this.kMvt; 
	}
	
	/**
	 * Обнуляет эелменты матрицы.
	 * @param A Матрица
	 */
	private static void zero(double[][] A) {
		for (int i = 0; i < A.length; i++)
			for (int j = 0; j < A[i].length; j++)
				A[i][j] = 0;
	}
	
	/**
	 * Разброс единиц в обобщённую матрицу.
	 * @param k Номер узла
	 * @param mvt Заданные перемещения
	 */
	private void throwOnes(int k, double[] mvt) {
		int i = 3 * k;
		if (bnd[0][i] > 0) // На диагонали перемещение уже учтено
			return;
		for (int j = 0; j < 3; j++, i++) {
			bnd[0][i] = 1;
			f[i] = mvt[j];
		}
	}
	
	/**
	 * Забрасывает перемещения в правую часть с соответствующими коэффициентами.
	 * @param r Номер узла
	 * @param Krs Матрица Krs
	 * @param sMvt Перемещения в узле s
	 */
	private void throwToRigth(int r, double[][] Krs, double[] sMvt) {
		int i = 3 * r;
		for (int k = 0; k < 3; k++, i++)
			for (int j = 0; j < 3; j++)
				f[i] -= Krs[k][j] * sMvt[j];
	}
	
	/**
	 * Забрасывает верхний треугольник матрицы Krs в обобщённую матрицу.
	 * @param r Номер узла
	 * @param s Номер узла
	 * @param Krs Матрица Krs
	 */
	private void throwTopTriangle(int r, int s, double[][] Krs) {
		int i = 3 * r;
		int j = 3 * s;
		for (int a = 0; a < 3; a++)
			for (int b = r == s ? a : 0; b < 3; b++)
				bnd[(j + b) - (i + a)][(i + a)] += Krs[a][b];
	}
	
	/**
	 * Уничтожает ссылки на ненужные экземпляры структур данных.
	 */
	private void clean() {
		crd = null;
		adj = null;
		mvt = null;
		D = null;
		coeff = null;
		Bi = null;
		Bj = null;
		Bm = null;
		Bp = null;
		B = null;
		Krs = null;
		kMvt = null;
	}
}